Stability of giant vortices in quantum liquids 
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We show how giant vortices can be stabilized in strong external potential Bose-Einstein con- 
densates. We illustrate the formation of these vortices thanks to the relaxation Ginzburg-Landau 
dynamics for two typical potentials in two spatial dimensions. The giant vortex stability is studied 
for the particular case of the rotating cylindrical hard wall. The minimization of the perturbed en- 
ergy is simplified into a one dimensional relaxation dynamics. The giant vortices can be stabilized 



only in a finite frequency range. Finally we obtain a curve for the minimum frequency needed to 

w-j ■ 

■ observe a giant vortex for a given nonlinearity. 
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: O Quantized vortices are characteristics of quantum systems. They appear for instance when a rotation 

is imposed to a quantum liquid. We study here the conditions under which vortices with multiple 

> ■ 

charges can be observed. It is in general not possible since the energy of an array of single vortices is 

o . 

^vq below the energy of a multiple vortex of same total charge. However an opposite trend appears for 
the energy if the system is confined. This in particular the case for the recent achieved Bose-Einstein 

o . 
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condensates where an optical trap maintains the gas in a small volume. We show that if the trap 
potential is strong enough the multiple charge vortex becomes stable. The linear stability of these 
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giant vortices is performed using one dimensional minimization techniques 

o 
o 



INTRODUCTION 
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has strongly renewed and restimulated the 



The recent achievement of Bose-Einstein condensates (BEC) 
numerous works and theories on quantum gases/liquids. By their lower densities, their larger healing lengthes, 
the BEC present important advantages compared to the traditional Helium superfluids, regardless the fundamental 
interest of the Bose-Einstein transition itself. One of the striking properties of these quantum gases/liquids is the 
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quantization of the vorticity ;4j]. It manifests through the formation of vortices whose circulation is quantified in units 
of h/m, where m is the mass of an atom of the gas and h is the Planck constant. Such vortices have been observed 



in BEC by two different procedures these last years|S|, |6|, |7| . Firstly 



ill 



by rotating the BEC through laser techniques. 



This experiment is the analogy of the famous rotating bucket of superfluid helium where the first direct observation 



of quantized vortices was obtained 



. The vortices appear suddenly as the rotation frequency increases, either one 



by one, or by group, depending on the experimental procedure. The other experiment (J 



7j exploits the dynamical 



transition to vortex shedding of an accelerated superflow 9]. Experimentally, to achieve the BE condensation the 
atoms of the condensate have to be placed in an optical trap which creates an effective (mostly harmonic) attractive 
potential. Thus, in the above experiments only single charged vortices were observed. No multicharged vortices 
(called giant vortex of charge q or q-vortex later on) were detected. However, by imposing specific initial states, it 
has been possible to observe metastable long-lived giant vortices experimentally ^j. It is in fact well known that 
multicharged vortices are always dynamically and thermodynamically unstable in homogenous and infinitely large 
quantum liquids n| . When the system is submitted to a trap potential, the situation becomes more tricky since the 
potential encourages vortex merging while the vortices try to expand the condensate. For BEC in harmonic trap, 
variational stability analysis suggests indeed the existence of a critical frequency above which the giant vortex should 
be energetically favoured to the single vortices lattice ^| . But this frequency is found higher than the frequency of 
the quadratic trapping potential. At such speed, the condensate is unstable since the centrifugal force overcomes the 
trap confinement. This explains why no stable giant vortices have been observed until now in rotating BEC, where 
only harmonic trap were experimentally considered. This peculiar property of harmonic confinement does not actually 
prevent the stability of giant vortices for stronger trap potentials. In fact different groups have shown numerically 
and argued analytically that giant vortices are stable for high enough rotation if the trap potential is taken stronger 
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than harmonic ll3l Il4l 1151 . Numerical simulations have exhibited these giant vortices and some coexistence diagramms 
were even obtained|l3l tl5]. although no clear stability arguments were presented. In the other work 14 the stability 
analysis was only performed in the framework of Wigner-Seitz approximation of a vortex lattice. Adequate strong 
trap potentials can be experimentally achieved thanks for example to Laguerre-Gaussian laser beams lfj but the 
observation of stable giant vortices remains still to be done. 

In this note, we address the stability analysis of the giant vortices for the rotating cylinder. We use numerical 
simulations for two-dimensional systems subject to rotation ft = Qz, assuming translational invariance in the z- 



3 

direction. The giant vortex formation is shown for illustration in a rotating BEC with quartic potential but we 
restrict thereafter our treatment to the case of an infinite cylindrical potential. 



THE MODEL 



The dynamics of a quantum Bose gas in the rotating frame in two dimensions is described through the dimensionless 



Gross-Pitaevskii (G-P) equation 



id t i> = f - - V 2 + V(r) + M\^\ 2 + iQdgj V (1) 

where t) is the condensate wave function normalized to unity: J dr\ip\ 2 = 1. r is the position vector (x, y), r = |r| 
and 9 the corresponding cylindrical coordinates. The nonlinear strength M = AirNa is related to the product of the 
particle number per unit length N with the s-wave scattering length a. This dimensionless equation has been in fact 
deduced from the usual model using a typical length scale d that is defined by the trap properties. Consequently the 
length, energy and time in equation Q are given in units of d, h 2 /(md 2 ) and md 2 /K (Ti the Planck constant and m the 
particle mass). The usual harmonic potential ^mLu 2 r 2 leads to V(r) = ^r 2 with d = y/K/JmuJt). Below we restrict 
our work to only two traps geometry. Firstly for illustration a quartic potential —jf—, so that V(r) = jr 4 and 
d = Wh/ (muj t ) again. This potential has a smooth behavior and describes a stronger trap potential than the usual 
harmonic potential of the experiments. The other potential models a quantum liquid entrapped in a cylindrical bucket 
of radius L. With d = L the potential is described by V(r) = if r < 1 and V(r) = oo otherwise. It corresponds 
to the well-known experiment on superfluid Helium where quantized vortices have been first vizualized 8] . Note that 
equation |JJ| admits only two control parameters for a given potential: the interaction coefficient M and the angular 
frequency f2. The dynamics derives from a free energy F through: idt'tp — Jzr with: 



F = Jdr QlWf + in^de^J + V(r) + yl^| 4 ) (2) 

Stable equilibrium solutions of Eq. correspond to local minima of F. The ground state is determined among 
the stable equilibrium solution as the one with the lowest energy . Stable equilibrium states of the condensate can be 
reached numerically through the dissipative dynamics of the Ginzburg-Landau equation. Formally it corresponds to 
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an imaginary time evolution of the GP equation (JTJ : 



dti/> 



( 




(3) 



We have introduced the chemical potential \i as the Lagrangian multiplier ensuring that the norm of ip is unity. 
Numerical simulations have been performed in two space dimensions up to 256 x 256 grid points. We have checked 
that numerical precision is not significantly affected by higher spatial discretizations for our range of parameters. 
Both Fast Fourier Transform (FFT) or finite difference methods have been used for spatial derivatives and again no 
dramatic changes were noticed between these schemes. The relative error in the normalization of the wave function 
was controlled to be below 10~ 3 and the precision for fj, is better than 10 -3 . 



To illustrate the formation of giant vortices, we follow numerically an equilibrium solution as CI evolves. We start at 
CI = with the ground state (vortex free solution) and we increase suddenly CI of a small increment SCI. We retrieve 
an equilibrium state through the imaginary time dynamics and the procedure is then successively reproduced. 
Both the quartic and the hard-wall potentials are considered. 

We first present the results for the quartic potential with M = 12tt (figures 0). Figure Q]a) shows the ground 
state solution, vortex-free at Cl = 0. This solution survives without any change until Cl = 5.4 where suddenly four 
vortices enter into the cloud (figure^b). They are well separated and form a square lattice. We observe also a spatial 
expansion of the cloud due to centrifugal effects. As CI increases further on(figure ^ c)), the four vortices approach 
slowly the center of the trap. In this case the full merging of the four vortices is attained at CI = 7.2 (figure ffld)), 
just before new vortices enter the cloud (at CI = 7.35). 

A similar situation is observed for a superfluid in a cylindrical trap although in that case the condensate cannot 
expand when vortices are nucleated. Indeed we account for the infinite hard wall potential by imposing as boundary 
conditions ip — for r > 1. Figure J5J) shows the density |V>| 2 of the solution of equation © as CI increases for 
M = 300. The nucleation of four vortices occurs at CI = 14.8 (figure (0 a)) from the vortex free solution. Further 
increments of CI lead to the convergence of the four vortices towards the center of the trap(figure J2|b)). But a 
different scenario arises here compared to the dynamics observed for the quartic potential. Indeed, before the four 
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Figure 1: Quartic trap potential: density plots \ip\ 2 for N ■ a = 3 at fl = a) 0., b) 5.4, c) 6.6, and d) 7.2. The grid size is 
128 x 128 and a pseudo-spectral method is used. The darker the color, the higher the density so that the vortices appear as 

white hole in the condensate. 

vortices collapse in a single 4- vortex, four new vortices are nucleated in the cloud at 17 = 20.8(figure (0 c)). The 
condensate is now composed of a lattice of eight well separated vortices. Then new increases of 51 lead eventually to 
the merging into a giant 8-vortex at Jl = 26.2 (figure (0d)). Above this rotation, when new vortices are nucleated, 
they are immcdiatly absorbed by the giant vortex. This persistence of giant vortices for increasing frequencies has 
been in fact observed for all the numerical simulations that we have performed. 
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Figure 2: Density plots for M = 300, for fi = a) 14.8, b) 18., c) 20.8, and d) 26.2. The grid size is 250 x 250 and the 
increment of the angular frequency is 0.2. The color scale goes from blue for zero density to red for the maximum density. 

The vortices appear as blue holes on the density field. 

STABILITY ANALYSIS 



By sake of simplicity, the stability analysis is performed for the cylindrical wall trap potential. We want to address 
the conditions under which conditions a q- vortex is a local minimum of the free energy @ . We determine firstly the 
giant q- vortex state by seeking solutions of the form: 
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^°(r) = f q {r)e^ 

The trivial symetry of the problem q — > — g with r2 > — makes us only consider later on the case q > 0. From 
the free-energy principle, to be an equilibrium solution f q has to satisfy: 



\ (/? + 7 - - M 4 3 + M M > n ) + « ■ = o (4) 

fi q (M, f2) is introduced again as the Lagrangian multiplier for the norm constraint < V'glV'q >=< /g|/g >== 
27r J rdrfq — 1, where < -|- > stands for the usual scalar product. The external potential term is replaced by 
the boundary condition -0 = for r > 1 so that this potential term is everywhere zero in the equation. We observe 
from the preceeding equation that the giant vortex solution is independant of £1 which only intervenes in the 
coefficient n q (M, O) = fJ> q (M) — q-Sl. These solutions are thus equilibrium states for all frequencies f2 regardless their 
stability. The free energy of the vortex writes: 



F q = 2n I r ■ dr ( 1((/') 2 + ^P) q Slf q + M f ) = T q (M) + U q (M) -q-Sl (5) 



Where T q (M) = tt J r ■ dr((f' q ) 2 + ipzf q ) is the kinetic energy and U q (M) — Mir J r ■ drf q the nonlinear energy 
term. E q (M) = T q (M) + U q {M) is the total energy of the vortex solution. The q- vortex function f q is found as before 
through the dissipative dynamics of the axisymetric Ginzburg-Landau equation: 



a/« = \ (/*' + 7 - - M ft + /*S/« (8) 

so that the equilibrium solutions f q are local mimina of the free energy F q . Figure JSJa)) shows the density f q of the 
solution for M = 40ir and q = 0, 1, 10 as a function of the position r. As it can be emphasized from equation 
the bigger q, the larger is the core of the vortex, and the higher is the density level outside the vortex core. Indeed, 
equation Q for r — > implies that the solution / g behaves like r 9 near the center. Figure 10b)) represents also the 
g = 5 solution for M varying between 4-7T and 8007T. The vortex core is shrinking as M increases in agreement with 
the evolution of the condensate healing length (oc l/\/M), so does the wall boundary layer. 
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Figure 3: Density of the condensate / 9 2 for a) M = 40% and q = 0, 1...10 b) q = 5 and M/pi = 4, 20, 40, 80, 200 400 and 800; 
the higher the maximum of the curve, the smaller M and the bigger q. The profiles have been obtained as the stationnary 
solution of equation © and a regular grid of only 100 mesh points is needed for convergence. 

Figure shows the free energy of the giant vortices at — as a function of q for different values of M. The 
energy is found convex with q for each curve. From this curve we can determine for each f2 which q- vortex minimizes 
the free energy. It defines an increasing set of critical frequencies fi™(M) such that the q- vortex is the one with the 
lowest energy for Q^_ 1 (M) < ft < ft™(M) with 



ft™(M) = E q+1 - E q 

A q- vortex is stable when it is a local minimum of the free energy J2J|. We expect from the numerics that such a 
vortex is unstable at low enough frequency ft: it decomposes itself into q single vortices (see figures Q and ©). On 
the other hand, high frequencies destabilize the g-vortex by nucleating new vortices which gather into a new giant 
vortex as observed for M = 100. Thus, the g-vortex can only be stable in a finite range of frequency and depending 
on the vaues of q and M, we may observe that the q- vortex can never be stabilized (it is the case for q — 4, M — 100). 
The stability analysis is performed by an expansion of the free energy around the g-vortex solution. We write the 
perturbed wave-function of the giant vortex for q > 1 in the form: 

Following the periodic dependance in 9, we decompose the perturbation 5ip q among azimuthal modes as a Fourier 
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Figure 4: Energy of the vortex solutions as a function of the charge q for three different values of M. From bottom to top, 

M = 1 (circles), M = 40 (squares) and M = 100 (diamonds). 
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The complex functions g^,(r) describe the q' pertubation mode of the q vortex and arc normalized to unity. The 
amplitudes e q > of the q' perturbation modes are complex numbers. We can omitt the q perturbation mode which will 
bring nothing to the stability analysis since ip® is already the minimizer of the axisymmetrical free energy with charge 
q. By analysing the perturbed wavefunction near r = 0, we observe that for q' < q the perturbation mode describes 
the destabilization of the q vortex into a g'-vortex in the center surrounded by q — q' single vortices. On the other 
hand, for q' > q the perturbation describes the nucleation from the wall of q' — q isolated vortices. We obtain for the 
perturbed wavefunction: 

i> = , = — 

The free energy correction SF writes at the first non-zero order in the e's: 
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sf= 530 





q'<q" 



\e q „J 



omitting the q superscript from now on. The set of operators £ q > 



,» follows: 



5F a , 



2irMfrf q (r) 2 g q ,(r)g q „(r)dr 



J* = 5 a > 



>q'+q",2q 



V 



2irM jrf q (r) 2 g* q ,(r)g q „(r)dr 



5F a „ 



J 



with the diagonal terms: 




(7) 



We notice here that the 9 dependance of the free energy has been integrated so that further calculations will only 
concern the r spatial variations. This axisymetrical simplification leads eventually to a one dimensional stability 
problem. The stability of the g-giant follows from the perturbated energy structure: if the operators C q > , q >> are a ^ 
found positive then the vortex is stable. On the other hand, if at least one of these operators is shown negative the 
giant vortex is unstable towards either decomposition into simpler vortices or nuclcation/annihilation of vortices. The 
operators £ q >, q " couple in fact the two modes g q / and gi q - q i so that the determination of their full spectrum remains 
a complicate task. 

However, if we could neglect the non-diagonal terms, the determination of the giant vortex stability would be 
straightforward. Indeed, the stability would be only determined by the minimization of the effective free energies 6F q / 
associated to the functions g q i for each q' ^ q. Following this assumption discussed in more details below, we only 
have to determine the eigenfunction g q i which minimizes the corresponding energy: 




(8) 



with the following conditions (omitting the hat of the functions g q > later on) : 



< g q ,\g q , >= 1. ^(0) = ,g^o(0) = g q ,{\) = 



(9) 
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We observe that energy J5J is unchanged through the symetry q' — ► — g' so that we just need to calculate the 
minimizcr for q' > 0. 5 9 ' corresponds in fact to the calculation of a q' giant vortex placed into the effective potential: 



U eff {r)=2Mf q {rf +V(r) 
The solution can be reached again through a dissipative dynamics: 

dtg g > = \ + 9 -f- jj<v) - 2Mf q g q , + ^g q , (10) 

where // is introduced here as the Lagrangian multiplier for the unity norm constraint. Figure JSJ shows precisely 
different g q i for q = 4 and M = 300. The q giant vortex profile is also indicated. As suggested by the effective 
potential, the perturbations modes are concentrated near the center of the trap for q' < q where the giant vortex 
density is very small. For q' > q the density g q > is concentrated near the wall since the centrifugal term of the 
perturbation dominates the g-vortex influence. Moreover, this figure is in agreement with our assumption of weak 
off-diagonal terms made above. Indeed, for |g'| < q the interaction term g(q')g(2q — q') is very small since one of the 
solution is concentrated near r ~ while the other one is only present near r = 1. However, for q' < — q the situation 
cannot be concluded easily since both term is concentrated near r = 1. 

By minimizing the energy SE q ^ q i (JHJ), we can determine whether the functionnal SF q > (JJJ is positive definite for all 
type of perturbations as function of the frequency O. The sign of this free energy for the g'-perturbation changes for: 



T q + 2U q - 5E q q * 



Q-Q 

It tells us that for q' < q the g' perturbation mode is stable for f2 > f2 c (g') while for g' > g it is stable only for 
fl < ri c (g'). The change of inequality is due to the change of sign of g — q'. We define then two critical curves 
Qiow{q') = ^c(g') and Qhigh(q') — ^ c (2g — g') calculated both for g' < g. For a given fi, if one can find a q' such 
that ^Iiow(q') > then the q- vortex cannot exist since it decomposes into more complex vortices configuration. On 
the other hand if one can find a q' such that SlhighW) < ^ the q vortex is also unstable and new vortices enter 
in the cloud. Thus the g-vortex can only be observed if Max g / <g f2; OUI (g') < Min 9 ' < q Qhigh(q') and in the range 
Ma,x q r <q fti ow (q') < il < Min q > <q n hi gh(q') the g-vortex is linearly stable. Figures JHJ a) and b) show £li ow (q') and 
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Figure 5: Density profiles of the perturbations modes g q i for a q — 4 giant vortex at M — 300. Only the profiles for q' = 
(dotted line), q' =2 (dashed line), q' = 6 (long dashed) and q = 20 (dot-dashed) are indicated. The q — 4 vortex profiles is 
also shown (solid line). The density for q' = and q' = 2 have been divided by height and two respectively for simplicity. 
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Figure 6: flhigh(q') (black circles) and Qiow (squares) as functions of the mode number q for M — 300, a) q — 4, and b) 
q = 8. It illustrates the two generic situations encountered for the stability of a (/-giant vortex: a) the 4-giant vortex is never 

stable; b) the 8-giant vortex is stable for Q in the range [25.09, 28.76]. 



Qhigh(q') for M = 300, q — 4 and q — 8 respectively. This figure is found in relatively good agreement with the 
full numerical simulations of the dynamics shown figure l|2*Jl. We remark on figure © that the 4- vortex cannot be 
stabilized since the perturbation mode q' — is still unstable when the high vortices perturbation q' = 18 becomes 
unstable (at ft — 22.29) as observed on figure (0). The solution with four single vortices becomes unstable at il = 
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Figure 7: Evolution of M&x q i <q Sli om (q') (circles) and Min g / <q Q.high(q > ) (red squares) for the giant vortex charge q between 2 
and 15. The two curves cross between q — 6 and q — 7, so that only above Q = 25.74 can giant vortices be observed. 
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Figure 8: Critical curves Q g (M) (black line) and q g (M) (red line). Giant vortices are linearly stable above these curves. 

which is close to the value found for the 4- vortex. For q = 8 

Figure (0 presents the evolution of Max g / <q Sli ow (q') < fi < Min q > <q Qhigh{Q') f° r Q S [2, 15]. As observed above, 
it shows that for q < 7 no giant vortex can be observed while for £1 > 25.74 giant vortices of charge q > 6 can be 
observed in a given range of f2. 

Defining Q g (M) (respectively q g (M)) as the lowest frequency (vortex charge) at which a giant vortex can be observed 
for a strength M, we determine in the f2-M plane (or q-M plane similarly) the region where the solution consist in 
stable giant vortices only. The curves Q, g (M) and q g (M) delimiting these regions are shown on figure JSJ). 
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DISCUSSION 

We have shown numerically that giant vortex solutions can exist in axisymmetrical 2-D BEC model for trap 
potentials steeper than harmonic. We have investigated in more details the linear stability of these vortices for the 
cylindrical wall-like potential. The axisymmetry property of the q- vortices simplifies the stability procedure into a one 
dimensional minimization problem. Under a no off-diagonal terms assumption presenting good agreement with the full 
numerical studies we can determine the regions where the giant vortices are stable. Finally a master curve presenting 
the frequency and vortex charge above which giant vortices are stable for varying nonlinear strength M is presented. 
These results have been obtained for the particular case of cylindrical wall potential but can be adapted to smooth 
BEC potentials. The stability analysis lead to the same minimization problem but the numerical determination of 
the q' modes needs more careful investigations since they are not concentrated on r < 1 disk. 

ACKNOWLEDGEMENTS 

It is my pleasure to thank Vincent Hakim for his encouraging remarks on this work. I also would like thank warmly 
the editors for organising this special issue. I cannot also finish without express my special thanks to Yves Pomeau 
who made me discover the magic world of nonlinear dynamics! 



[1] C.C. Bradley, C.A. Sackett, and R.G. Hulet, Phys. Rev. Lett. 75, 3969 (1995). 

[2] M.H. Anderson, J.R. Ensher, M.R. Matthews, C.E. Wieman, and E.A. Cornell, Science 269, 198 (1995). 
[3] K.B. Davis, M.-O. Mewes, M.R. Andrews, N.J. van Druten, D.S. Durfee, D.M. Kurn, and W. Ketterle, Phys. Rev. Lett. 
75, 3969 (1995). 

[4] R.J. Donnelly, Quantized Vortices in Helium II (Cambridge University Press, Cambridge, England., 1991). 
[5] K.W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Phys. Rev. Lett. 84, 806 (2000). 

[6] C. Raman, M. Khl, R. Onofrio, D. S. Durfee, C. E. Kuklewicz, Z. Hadzibabic and W. Ketterle, Phys. Rev. Lett. 83, 2502 
(1999). 

[7] R. Onofrio, C. Raman, J. M. Vogels, J. R. Abo-Shaeer,A. P. Chikkatur and W. Ketterle, Phys. Rev. Lett. 85, 2228 (2000). 
[8] E.J. Yarmchuk, M.J.V. Gordon, and R.E. Packard, Phys. Rev. Lett. 43, 214 (1979). 
[9] T. Frisch, Y. Pomeau and S. Rica, Phys. Rev. Lett. 69, 1644 (1992). 
[10] P. Engels, I. Coddington, PC Haljan, V. Schweikhard, and E.A. Cornell, cond- mat/0301532 (2003). 



15 

[11] I.S. Aronson and V. Steinberg, Phys. Rev. B 53, 75 (1995). 

[12] Y. Castin and R. Dum, Eur. Phys. J. D 7, 399 (1999). 

[13] E. Lundh, Phys. Rev. A 65, 043604 (2002). 

[14] U.R. Fischer and G. Baym, cond-mat/0111443 (2001). 

[15] K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. A 66, 053606 (2002). 

[16] T. Kuga, Y. Torii, N. Shiokawa, T. Hirano, Y. Shimizu, and H. Sasada, Phys. Rev. Lett. 78, 4713 (1997). 

[17] E.P. Gross, J. Math. Phys. 4, 195 (1963). 

[18] L.P Pitaevskii, Sov. Phys. JETP 13, 451 (1961). 

[19] F. Dalfovo, S. Giorgini, L. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999). 



